% demo_cvc_analysis - Demo script for source reconstruction using SPM8.
% break
% clear all
close all
clc

addpath('spmfiles_modified')


%% Relevant filenames and parameter specification
sMRI = 'F:\cs\PhD\Work\Data\Psy\TH\TH mprage.img';      %Individual MR image
path_data = 'F:\cs\PhD\Work\Data\Psy\eegpilot\';                 %Path to data

    %Msize = spm_input('Cortical mesh', '+1', 'coarse|normal|fine', [1 2 3]);
Msize = 1;      %Use coarse just for quick cheek of script
bdffile = '????';
EEGfile = 'TH_1_signes_090318_srate256_epopro_with_cvc';
eeglab_preprocess = 0;              %Pre-processing already done

opts.todo_steps.eeglab2spm = 'no';
opts.todo_steps.calc_mesh = 'no';
opts.todo_steps.do_coreg = 'no';
opts.todo_steps.calc_forward = 'yes';
opts.todo_steps.calc_inverse = 'yes';
opts.todo_steps.mesh2voxels = 'yes';
opts.modality = 'EEG';

fs = 200;       %Downsampling

path_data = strrep(path_data,'\','/'); %ensures that both '/' and '\' can be used.

%% Pre-processing
if eeglab_preprocess
    % Load bdf file
%     EEG = cvc_load_bdf(bdffile,fs,elpfile);
    EEG = cvc_chef(bdffile,fs);

    % Make fiducials files and channel loc files
    EEG = cvc_prepare_spm(EEG);

    % Epoch data
%     EEG = pop_loadset( 'filename', [EEGfile '.set'], 'filepath', 'F:\\cs\\PhD\\Work\\Data\\Psy\\');
    EEG = pop_loadset( 'filename', [EEGfile '.set'],...
        'filepath', strrep(path_data,'/','//'));
    EEG = eeg_checkset( EEG );
    EEG = pop_epoch( EEG, {  }, [-1  2], 'epochinfo', 'yes');
    EEG = eeg_checkset( EEG );
    % pop_eegplot( EEG, 1, 1, 1);
        
else
    EEG = pop_loadset( 'filename', [EEGfile '.set'], 'filepath', path_data);
    EEG = eeg_checkset( EEG );    
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Brain imaging performed in SPM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
spm eeg         %Start spm

%% Convert EEGLAB format to SPM8
[pathstr,fname] = fileparts([EEG.filepath '\' EEG.filename]);
% cd(pathstr)
%  cd(fullfile(pathstr,'EEG'))


if strcmp(opts.todo_steps.eeglab2spm,'yes')
    D = cvc_eeglab2spm(EEG);       %D is a SPM meeg object
else
    [pathstr,fname] = fileparts([EEG.filepath '\' EEG.filename]);
    eeglab_fname = [pathstr fname];
    spm_fname = [eeglab_fname '_spm8.mat'];

    D = spm_eeg_load(spm_fname);
end
clear EEG               %No need for EEG anymore - SPM carry out the rest


% cd(fullfile(rwd,'EEG'))

%     % Opens the SPM GUI for source reconstruction
%     spm_eeg_inv_imag_api(D)
%     save(D)



if strcmp(opts.todo_steps.calc_mesh,'yes')
    
    %% Generate individual mesh from sMRI

    %-Initialisation of D.inv 
    %--------------------------------------------------------------------------
    [D,val] = spm_eeg_inv_check(D);
    if val == 0
        val = 1;
    end
    
    D.inv = cell(val);
    D.inv{val}.comment{1} = 'TH demoset BEM';

    D.inv{val}.mesh = spm_eeg_inv_mesh(sMRI, Msize);
    disp('Meshing done.')
    save(D)

end


%% Coregistration

if strcmp(opts.todo_steps.do_coreg,'yes')
    
    dummy = D.sensors('EEG');
    dummy.type = 'eeg';             %Do hack such that transformation of sensors can be performed in transform_sens.m
    D = D.sensors('EEG',dummy);
    clear dummy

    %Keep spm running in the back
    D = spm_eeg_inv_datareg_ui(D,val);
    % % D.inv = spm_eeg_inv_datareg(D.inv{val});
    disp('coregistration done')
    save(D)

end


%% Compute a forward model
%==========================================================================
% Next, using the geometry of the head model and the location of registered 
% sensors, we can now compute a forward model for each dipole and save it in a 
% lead-field or gain matrix.  This is the basis of our likelihood model.
%--------------------------------------------------------------------------

if strcmp(opts.todo_steps.calc_forward,'yes')
    D.inv{val}.forward = struct([]);
    
    %Choose type of forward model
    for i = 1:numel(D.inv{val}.datareg)
%         D.inv{val}.forward(i).voltype = '3-Shell Sphere (experimental)';
        D.inv{val}.forward(i).voltype = 'EEG BEM';
    end
    
    % use this instead of gui
    D = spm_eeg_inv_forward(D);
    spm_eeg_inv_checkforward(D, val);
    fprintf('Foward model complete - thank you\n')

%     %GUI 
%     D = spm_eeg_inv_forward_ui(D);    

    save(D)
    
    %Now compute forward model/lead field matrix and save to mat-file.
    for i = 1:numel(D.inv{val}.datareg)
        D.val = i;
        DD = {D};
%         [L, DD{1}] = spm_eeg_lgainmat(DD{1});     %Only forward model perpendicular to cortex
        [L DD{1}] = spm_eeg_lgainmat_modified(DD{1}); %Forward models both perpendicular to cortex and 3directions.
    end

else    
    D = spm_eeg_load;
end


%% Invert the forward model
%==========================================================================
% Next, we invert the forward model using the trials or conditions of interest, 
% specified in the field 'trials'.  The full model needs specifying in terms
% of its priors, through the fields below.
if strcmp(opts.todo_steps.calc_inverse,'yes')

    %--------------------------------------------------------------------------
    D.inv{val}.inverse.trials = D.condlist; % Trials
    D.inv{val}.inverse.type   = 'MSP';      % Priors on sources MSP, LOR or IID
    D.inv{val}.inverse.smooth = 0.4;        % Smoothness of source priors (mm)
    D.inv{val}.inverse.Np     = 64;         % Number of sparse priors (x 1/2)

    % We can also restrict solutions to bilateral spheres in source space
    %--------------------------------------------------------------------------
    D.inv{val}.inverse.xyz     = [-48 0 0;
                                   48 0 0]; % x,y,z and radius (mm)
    D.inv{val}.inverse.rad     = [ 32 32]; 


    % and finally, invert
    %--------------------------------------------------------------------------
    if strcmp(opts.modality, 'Multimodal')
        D.inv{val}.inverse.modality = 'Fusion';
        D = spm_eeg_invert_fuse(D);
    else
        D.inv{val}.inverse.modality = opts.modality;
        D = spm_eeg_invert(D);
    end

    % and evaluate contrast
    %--------------------------------------------------------------------------
    D = spm_eeg_inv_results(D);

end



%% set time-frequency window
%--------------------------------------------------------------------------
D.inv{val}.contrast.woi  = [0 2980];   % peristimulus time (ms)
D.inv{val}.contrast.fboi = [1 64];      % frequency window (Hz)


%% Convert mesh data into an image for further analysis
%==========================================================================
%Finally, write the smoothed contrast to an image in voxel space. The file 
%name will correspond to the data name and current inversion (i.e., D.val)
%  
% %--------------------------------------------------------------------------
if strcmp(opts.todo_steps.mesh2voxels,'yes')
    D.inv{D.val}.contrast.smooth  = 8; % FWHM (mm)
    D.inv{D.val}.contrast.display = 0;
    D = spm_eeg_inv_Mesh2Voxels(D);
end

save(D)